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B I Abstract 

^ ■ This paper reconsiders the problem of calculating the expected set of probabilities (pi), 

^ ■ 

given the observed set of items {rrii}, that are distributed among n bins with an (unknown) 

o : 

1^ ■ set of probabilities {pi} for being placed in the ith bin. The problem is often formulated 

o ' 

' using Bayes theorem and the multinomial distribution, along with a constant prior for the 

^ , values of the pi, leading to a Dirichlet distribution for the {pi}. The moments of the pi 

'X: 

^ ' can then be calculated exactly. Here a new approach is suggested for the calculation of the 

moments, that uses a change of variables that reduces the problem to an integration over 
a portion of the surface of an n-dimensional sphere. This greatly simplifies the calculation 
by allowing a straightforward integration over (n — 1) independent variables, with the 
constraints on the set of pi being automatically satisfied. For the Dirichlet and similar 
distributions the problem simplifies even further, with the resulting integrals subsequently 
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factorising, allowing their easy evaluation in terms of Beta functions. A proof by induction 
confirms existing calculations for the moments. The advantage of the approach presented 
here is that the methods and results apply with minimum or no modifications to numerical 
calculations that involve more complicated distributions or non-constant prior distributions, 
for which cases the numerical calculations will be greatly simplified. 

1 Introduction 

Many problems involve placing N objects into n bins, with probabilities pi for the object being 
placed into the ith bin. Given the values of the set of {pi}, the probability density P(mi, m2, m„|pi,p2, •••,Pn) 
for the distribution of the set of {rrii} objects can be calculated, and is well-know as the multi- 
nomial distribution, 

P(mi,m2,...,m„|pi,p2,--,Pn) = j j ;'^7=iPT" (1) 

with the constraint that ELiPi = 1 and T,7=im = N. Bayes theorem, P{A\B)P{B) = 
P{B\A)P{A) requires, 

P{pi,P2, ...,p„|mi,m2, ...,m„)P(mi,m2, ...,m„) = P(mi,m2, m„|pi,p2, ...,Pn)P{Pi,P2, ■■■,Pn) 

(2) 

that in principle allows us to calculate P(pi,p2, ^2, •••"^n). the probability of the set 

of probabilities {pi} with i = 1 to i = n, given the observed set of {rrii}. Often in such prob- 
lems, P{pi,p2, ...,Pn) is taken to be constant, and P(mi,m2, ...,m„) is chosen to ensure that 
P{pi,P2, ...,Pn) is correctly normalised [T]. Applying this approach to the multinomial distribu- 
tion, leads to a Dirichlet distribution, for which exactly calculated moments can be obtained. A 
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recent approach to this problem by Friedman [2], relied on an identity discovered by Gauss that 
involves the integral representation of the hypergeometric distribution. The same is true of a 
recent exact calculation that corrects conjectured but widely used mark-recapture estimates [3J, 
this and the coincidental timing of its revision on arXiv are what brought this problem to the 
author's attention. 

Here an alternative method of calculation is considered. I suggest a change of variables that 
elegantly leads to a simple calculation for the moments of the {pi}, and confirms existing results. 
The advantage of the method is that it can be applied very generally, and allows comparatively 
straightforward numerical integrations for the most general situations when analytical solutions 
may not be possible. The crux of the problem is the integration of a function over all possible 
values of Pi between and 1, subject to the constraint of YH=iPi = 1- This appears in many 
situations, the specific case considered here is the product n"^^p™' that arises in the Binomial, 
Multinomial, and Dirichlet distributions for example. 

2 The Calculation 

Consider the integration of the product 11"^^^™^% over all sets of values of the pi, subject to the 
constraints of < < 1 for all i, and Y.7=iPi = 1- Casella and Berger [4J, the moments 
are obtained by a delightful trick (page 181), that simplifies the problem to integration over a 
binomial distribution. In Friedman [2] the integral is accomplished by a nested set of integrals, 
each of which depends on the calculation of the integrals within it, with for n = 3 for example. 




(3) 
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where Y.'i=iPi = 1 has been used to write P3 = I —pi —p2- Here I start in a similar way, writing, 

^UpT = (i - E j ^'l=lvT (4) 

that for n = 3 is - P\- p-i)^'^ ■ Eq. H recognises that the constraint of YIi=\Pi = 1 

leads to (n — 1) free parameters, or 2 free parameters for n = 3. For a radius of r = 1 the 
n-dimensional polar co-ordinates are: 

Xi{n) = cos 01 

X2{n) = sin 6*1 cos 6*2 

X3 (n) = sin 61 sin 62 cos 6'3 

(5) 

Xn-i{n) = sin ^1 sin 6'2... sin 6'n_2 cos 
Xn{n) = sin6'i sin6'2... sin6'„_2sin6'„_i 
Notice that Xi{n) and Xi{ny will vary continuously between and 1 as the set of 6i are varied 
continuously between and n/2. Also notice that SILi^il"^)^ = 1' ^""^ consequently that 
Xnin)"^ = 1 — Y^i=lxi{ny. Therefore the substitutions of pi = xi{n)'^, p2 = X2{nY, ... , 
Pn-i = Xn-i{ny, will ensure that Yli=iPi = 1. integrals over 9i from = to n/2 will 
allow Pi to vary continuously over all values between and 1. 

Note that the constraint of Z^ILiPi = ^ \eads to (n — 1) free parameters, that after the 
change of variables, correspond to the set of 9i with i = Ito (n — l). Also note that although we 
are using polar co-ordinates in n dimensions, because we have set r = 1, there are only (n — 1) 
free parameters. 

The Jacobian of the co-ordinate transformation is J = \dxi{ny /d9j\. Notice from Eq. [5] 
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that dxi{nY /dOj = for j > i. Consequently the determinant has zeros above the diagonal, 
and will evaluate easily to give J = n"J"/ \dxi{ny /d9i\. 

Before proceeding to the general case, consider again the case with n = 3, for which case. 



xi(3) = cos 6*1 
X2(3) = sin 6*1 cos 6*2 
0:3(3) = sin 6*1 sin 62 
The product (l — Pi) " ^7=iPT^ becomes, after the change of variables, 

(1 - Pi - P2r' pTpT = (sin' Oi sin' ^2)"' (cos' ^1)™^ (sin' 9^ cos^ 62 

= (cos^'"! 01 sin'("^2+'"^) ^1) (cos2™2 0^ sin""3 a 

The Jacobian is. 



(6) 



7712 



(7) 



J 



(8) 



—2 cos 61 sin 61 
2 sin 9i cos ^1 cos^ 62 —2 sin' ^1 sin 62 cos ^^2 
1^2 cos 61 sin^ 6'i) (2 sin 62 cos 612) 
Therefore using Eqs. [7] and [5] the integral in Eq[3]can be equivalently calculated from. 



/3 = /"^^ dOi j^'^ dd2 (cos""^ Oi sin'("2+™3) (^cos'™2 ^2 sin'™« ^2) (2 cos Oi sin^ ^1) (2 sin 62 cos ^2 







(9) 



This integral factorises into. 



/3= 1^2 |^"^'rfeiCos'(™^+i)-i0isin'("^^+'"^+')-i^^i^ [2 |^"^'de2cos'('"^+i)-i02sin'("^^+i)-i^^2) 

(10) 

the above Eq. (TUlwill be used as a starting point for a proof by induction for the general case 



later. 



Many readers will immediately recognise the integrals as Beta functions, and it is well known 
that, 

2 r^'d^cos^^-i^sin^"-!^ = B(m,n) = gl^^lE^ (n) 
Jo ^ ' r{m + n) ^ ' 

Consequently /a is easily evaluated as, 

^ r(mi + l)r(m2 + mg + 2) T^m^ + l)r(m3 + 1) 
r(mi + m2 + ms + 3) r(m2 + mg + 2) 

Cancelling terms and writing in terms of factorials this gives, 

h = , (13) 

(mi + m2 + ms + 2)! 

For non-integer Eq. [TT]must be left written in terms of Gamma functions. 

If we now wish to calculate (pi) for example, we simply need to evaluate /3(mi+l, m2, m^) / I^{mi, m2, m^) 
(mi + 1)/ (mi + m2 + + 3) = (mi + 1)/(A^ + 3) with = mi + m2 + m^, as found by 
Friedman. Other moments are easily calculated in a similar way. 

For the general case, consider the formulae, 

/■7r/2 1-71/2 i-k/2 

In = dOij^ de2... den-iU]llKj{n) (14) 
Kj{n) = 2 cos2("^+i)-i(^^,) sin'^"=.+i('+"')~'(0,) (15) 

where I note that I]r=i+i(l + "^«) = ('^ ~ j) + I]r=j+i"^^ ^"'^ the dependency on n of Kj{n) is 
through the upper limit in the sum. Note that Eqs. [T4]and [T5]are true for = 3, as can be seen 
by comparison with Eq. [101 I will assume this is true for n = k then show that this implies it is 
true for n = + 1, and consequently for all A; > 3 by induction. 
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Firstly consider the integral with n — k. For n — k the change of variables is, 

Pi — Xi{ny =cos^^i 

P2 — X2{ny = sin^ 9i cos^ 6*2 

P3 = (n)^ = sin^ 9i sin^ 62 cos^ ^3 

= Xk-i{ny = sin^ 6'i sin^ 62... sin^ 6'fe_2 cos^ 9k-i 
Pk ^ Xk{nY = sin^ 6*1 sin^ 6*2 .. . sin^ 6'fc_2 sin^ 9k-\ 
and the integrand is 11^^ ^p^*, with a Jacobian that as noted previously, simplifies 
li^iZl \d{xi{kf)/d9i\. This gives the integral 4 as, 

h = j^'^ dOi... j^'^ dOk-iTltMkf'^^Ii'^zl \dxj{kf/d9_ 

Now consider n — k-\-l, for which the change of variables is. 

Pi — Xi{nY —cos^Oi 
P2 = X2{ny = sin^ 9i cos^ 62 
Pz — Xz{nY — sin^ 9i sir? 62 cos^ ^3 



Pk-i = Xk-i{ny = sin^ 9i sii? 62... sii? 9k-2 cos^ 9k-i 

Pk — Xk{nY = sin^^isin^ ^2 ■■ -8111^^^-2 sin^^fe-icos^^fc 
Pk+i = = sin^ 9i sir? 92... sii? 9k-2 sir? 9k-i sir? 9k 



and the integral Ik+i is. 



4+1 = j^'^ d9i... d9kli\llx,{k + l)2™»n5=i \dx,{k + lf/d9^ 



Now notice that for i = 1 to i = (A; — 1), Xi{k) = Xi{k + 1). For i = k, Xk{k + 1) = Xk{k) cos^ 9k- 
Therefore, 



Wl^^Xiik + = ntiX,(A;)2-« cos2-'=(^fc)xfe+i(A; + 1) 



,2m 



k + l I 



Similarly the Jacobian can be written as, 



Therefore we have. 



rn/2 



■k/2 



\2TTfc-l 



sin2('"*+i+i)(^^)...sin2(™*+i+i)(^fe„i)2cos2(™''-+i+i)-i(4)sin2(™^+i+i)-i(efc) 



Comparing Eq. [IT] with the assumption of Eq. [141 we find. 



dxi{kf 



(20) 



(21) 



(22) 



(23) 



with Kj{k) given by Eq. [151 Under this assumption the integrand of Eq. [221 can be written as, 



2 cos2('"'=+i+i)-i(0fc) sin2('"'=+i+i)-^(4^ 



2(mfc+i + l) 



(24) 



Note that. 



Kj{k) sin2(™'=+i+i)(^j) 



The extra factor in Eq. [24l is. 



2cos2(-^+i)-i(%)sin'S'=^+i('+'"')-'(^^,-) 
+ 1) for 1 < j < (A; - 1) 



(25) 



2pQg2(mfe+i+i)-i(^^^) sin2('"^-+i+i)-i(4) = Kfe(A; + 1) 



(26) 



Therefore we have, 

h+i = de,... r'^ de^ut^K.ik + 1) (27) 

JO JO 

which is just Eq. [14] with n = (A; + 1), and Ki{k + 1) as given by Eq. [151 Since we've shown 
Eq. [27] to be true for n = 3 and that its truth for n = k implies it to be true for n = {k + 1), 
then by induction Eqs. [14] and [15] are true for all n>3. 

Eq. [2Z] is easy to evaluate. Because 9i only appears in Ki{k + 1), the integral factors into, 

4+1 = nil r'\d,K,{k + i) (28) 

JO 

Noting Eq. [T5]for Ki{k + 1), each of the integrals can be recognised as a Beta function, with. 



(29) 

where in the denominator of the last line we used + 1 + Z]f=j+i(l + "^/) = Z]f=/(1 + "^O- To 
obtain an explicit value for the integral, now we simply need to multiply out the terms, with. 



^ r(mfc_i + l)r(mfc+mfc+i+2) ^ r(mfc+l)r(mfc+i + l) 
r(r)ife_i+m,fe+rrtfc+i+3) r(mfc+mfc+i+2) 



(30) 



Cancelling successive terms, leaves. 



_ r(mi + l)r(m2 + l)...r(mfc + l)r(mfc+i + l) 



which when written in terms of factorials and = YaIi^^u gives 



For non-integral values of mj the Eq. [32] must be remain expressed in terms of Gamma functions. 
Note that the above expression (13^ is for n = A; + 1, and usually we will evaluate it with n = k, 
for which case Ik = mi\m2\...mk\/{N + k — 1)!. 

To obtain the gth moment of pi one simply needs to substitute (m^ + q) for nii in I^, and 
calculate the ratio of /^(mj + q)/ Ii^^nii), whose meaning is hopefully clear. For example, {pi) is 
given by, 

_ mi\m2\...{mi + l)!...mA;! {N + k - 1)! _ + 1 

~ {N + k)\ mi\m2\...mk\~ N + k ^ ' 

where the notation {pi) is used to denote the moment of when there are k "bins". Similarly, 

, 2x ^ mi\m2\...{m, + 2)!...mfc! {N + k - 1)! ^ (m^ + 2)(m^ + l) 

~ (iV + A; + l)! mi\m2\...mk\ ~ {N + k + l){N + k) ^ ' 

Giving the standard deviation as, 

^^^^ ^^'^ ~ (N + ky{N + k + l) ^ ' 

These results are in agreement with those of Friedman. Higher order moments are also easily 
calculated. The difference of the skewness from zero for example, can give an indication of 
the extent to which noise in the data should be regarded as non-Gaussian. Note that because 
J:7=lP^ = 1. then. 



^ ^ (36) 

= i:ii(p.) 

where D is used as shorthand to indicate that the integral should be over the correct domain of 
integration subject to the constraint of YJi=iPi = 1- Eq. |36] is correctly satisfied by Eq. 
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3 Remarks 

There are a variety of distributions in which the {pi} only appear in a factor of II^^^p^S and 
the results here apply to those cases also. More generally the probability distribution or its prior 
could involve any function of {pi}- For example, we might want to introduce a suitable prior 
into the problem so as to bias against "outliers", or towards a particular set of {pi}. In these 
more general cases the change of variables to n-dimensional spherical polars will still allow a 
comparatively straightforward numerical integral. A numerical integral over the {pi} subject to 
< Pi < 1 and J2i=iPi — 1, without the change of variables to spherical polars, is not so easy. 
For some combinations of priors and probability distributions the integral will remain factorisable 
after the change of variables. This might continue to be useful for other analytical calculations. 
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